Trace formula for counting nodal domains on the boundaries of chaotic 2D billiards 
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Given a Dirichlet eigenfunction of a 2D quantum billiard, the boundary domain count is the number 
of intersections of the nodal lines with the boundary. We study the integer sequence defined by these 
numbers, sorted according to the energies of the eigenfunctions. Based on a variant of Berry's random 
wave model, we derive a semi-classical trace formula for the sequence of boundary domain counts. The 
formula consists of a Weyl-like smooth part, and an oscillating part which depends on classical periodic 
orbits and their geometry. The predictions of this trace formula are supported by numerical data computed 
for the Africa billiard. 
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Recently, the study of nodal patterns witnessed a remark- 
able renaissance, and attracted the active interest of scien- 
tists from very diverse fields — quantum chaos, acoustics, 
optics, spectral theory, percolation and more [1]. The num- 
ber of nodal domains — connected components on which a 
real wave unction has a constant sign, is an important fea- 
ture, which has been used to characterize eigenfunctions 
of wave equations and even to resolve iso-spectral ambigu- 
ities [0]. For quantum billiards, the number jy^ of nodal 
domains of the n'th eigenfunction (sorted by increasing 
eigenvalues En < ^n+i) can never exceed n [3]. The 
number of nodal domains can be computed explicitly for 
separable systems. However, in the non-separable cases, 
there exists no analytical tool which provides the number 
of nodal domains, and even the numerical counting prob- 
lem is difficult due to the dependence on the detailed struc- 
ture of the domains. In pi] it was shown that the asymptotic 
distribution of depends on the dynamics of the underly- 
ing classical system. In the separable case, the parameters 
of the classical phase space determine the nodal counts in 
the semi-classical limit. In the chaotic case, the distribution 
matches the predictions of a percolation model, which was 
proposed (but not rigorously justified) in [5]. 

The partition to nodal domains induces a partition of the 
boundary to "boundary domains". They are defined as the 
nodal domains of a"boundary function", which for Dirich- 
let billiards, is the normal derivative of the eigenfunction 
on the boundary. In 2D, the number rjn of boundary do- 
mains equals the number of nodal points separating them — 
the boundary intersection (BI) points. This number is more 
accessible, both numerically and theoretically, than z/^. At 
the same time, it carries the fingerprints of the underlying 
classical dynamics of the billiard [4]. Also, the number 
of boundary domains provides information which is essen- 
tial for estimating the total number of nodal domains 
It was recently shown [7] that rjn = 0{y^). In [0], the 
asymptotic distribution of rjn/ y/n was computed for a class 
of integrable billiards. For a chaotic billiard of area A 
and boundary length £, a random wave model yields the 
estimate [S S] ^ >C^/(27r), where q = ^/^im/A is 



the leading asymptotic estimate for the n'th wave-number 
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The purpose of the present paper is to go beyond the 
simple estimate r]n ^ £g/(27r), and provide a trace for- 
mula which approximates the mean value as well as the 
fluctuations in the sequence T\n, in terms of the periodic or- 
bits of the classical billiard. Such formulae were proposed 
in the past for the counting of nodal domains in separable 
billiards [|9i]. Here is the first time that a counting trace 
formula is written down for the chaotic case. 

The counting of boundary intersections is performed by 
computing the density dj^{n) = XlmGN* ^(^ ~ T^)Vm' For 
a chaotic Dirichlet billiard with a smooth boundary, we de- 
rive in the sequel the following asymptotic expression: 
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where p enumerates classical periodic orbits, r e N* 
counts repetitions of the orbit, Lp is the length of the or- 
bit, Mp the monodromy matrix, Up the Maslov index, and 
q = q + C/ {2A). <I?p is a trigonometrical factor depending 
on the Up bounce angles of the orbit p: 
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Before sketching the derivation of the trace formula, we 
shall demonstrate its application. We computed the low- 
est 20,000 eigenfunctions (0 < /c^ < 260) of the Africa 
billiard [floll . the corresponding BI count sequence 77^, 
and the density df^{q) = p{q — ^/4:7m/A)r]n (where 
here and in what follows p{x) is a narrow Gaussian ap- 
proximating the Dirac 5. p is defined as a density, so 
p{q — Qo) = — ^0) • w4go/(27r)). Subtracting the 
predicted smooth part and scaling, we computed 

f{q) = (d^ — d^™) /q'W{q), where is a Gaussian "win- 
dow function" of width a = 50 and center qo = 130, which 
was used for softening the sharp cutoff due to the finiteness 
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of the computed spectrum. The length spectrum, which 
is the Fourier transform f{x) of f{q) was compared with 
fsc\{x)^ the theoretical prediction based on the oscillating 
part of ([T]). The latter was computed using 70 classical pe- 
riodic orbits (with up-to 7 bounce points) and 15 complex 
periodic orbits whose lengths had a very small imaginary 
part (however, the complex orbits did not have a significant 
effect). 

FIG. [Ha) displays several peaks centered at lengths of 
periodic orbits which match quite well with the theoretical 
predictions. A more detailed comparison is presented in 
FIG.[TJb). Due to (O, orbits whose angles are close to 60° 
are inhibited. Indeed, the triangular periodic orbits of the 
billiard, whose lengths are in the range 5.07-6.05, cannot 
be seen above the background level. The structure around 
X = 6.5 is due to several periodic orbits that pass very close 
to the boundary at the region of its highest concavity. The 
poor agreement between the semi-classical theory and the 
numerics in this region is due to penumbra corrections [[TTll 
which were not included. The random background of am- 
plitude ^0.04 observed in the plots does not seem to di- 
minish when qo increases (within the range of our study). 
This phenomenon will be discussed below. 

The derivation of the trace formula ([T]) starts by express- 
ing r]n Sisrjn = § bn{s)ds, where < 5 < £ is the bound- 
ary arc length and the BI density bn{s) is given by 
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5{Un{s)) \Un{s)\. (3) 



and s^""^ are the BI points 



Here, Un{s) = n{s) • V'0n(^(-§))/^n is the scaled normal 
derivative of the n'th eigenf unction t/jn, taken at the bound- 
ary point r(s), Un{s] 
(zeros ofun)' 

The trace formula will be derived for a smoothed version 
d^{n) J drj{m)p{n - m) dm = EmGN* rn)rjm, 
(The smoothing kernel p was defined above). In the sequel 
we shall consistently use the notation {y) to denote the 
p smoothed density of the quantity X in the variable y. 

Consider the spectral density of 6 (s) at wavenumber k, 
di){s]k) = ^n5{kn — k)bn{s). It is approximated by 
<(s;fc) = dP{k){bn{s))n, where d^(fc) = En P(^n - k) 
is the smoothed spectral density and {b{s))n = En ^n{s) • 
p{kn — k)/d^{k) is the spectral average of bn{s) around 
k. If we choose the width of p to be of order /c~2, then 
as ^ oo, the discrete ensemble of boundary functions in 
the corresponding spectral window around k approaches a 
continuous distribution. To proceed, we introduce the con- 
jecture that for chaotic billiards this limiting distribution 
is Gaussian. This conjecture, may be seen as a variant of 
Berry's random wave model JT^l, adjusted to the boundary 
as in I13I] and (l^ . While the validity of the conjecture is 
expected to improve in the semiclassical limit, numerical 
tests within the range of k values used here reveal a resid- 
ual Kurtosis which might explain the background observed 



in the length spectrum. Adopting the random waves con- 
jecture enables us to express the mean of the density b{s) 
in terms of the variances of the field u{s) and its derivative 
iiis) 111: 
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To compute the required densities, we write 

E2k 
5{kn - k)un^{s) — Im5f(s, s; k) (5) 
TT 

n 

where g{s,s']k) = Y^j^Un{s)un{s')/ [kJ - /c^) is the 
boundary Green function. As shown in JT^, g can be ex- 
panded as ^ = ^^9^^ where h is the integral oper- 
ator with kernel function h. The functions qq and h are 
given by 

h{s, s'] k) = 2n(s) . V^(,)Go (r(5), r(sO; k) 

2^ , d^G^{r{s),r{s')',k) 
go{s,s ]k) = —2^n,n j- 
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and Go = \H'^{k\r — r'\) is the free Green function in 
2D. To handle the singularities involved in this expansion 
and make it more amenable to semi-classical treatment, we 
choose a large cutoff 1 <C <C kC, and split g^) into a 
"far" (off diagonal) part and a "near" (close to the diagonal) 
part: 

gi^^ = 9q{s, s')H{xc/k -\s- s'\), 
9o^^ = 9o{s, s)H{\s - s\ - xc/k) 



(where H is the Heaviside step function). It can be shown 
that for large k 
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where f denotes Hadamard finite part integration the 

near diagonal part g[^'^ is bounded and hg[^'^ is negligible. 
Hence, the expansion of g can be rewritten as: 



n=0 



(6) 



Substituting this result in (|5]), we get a similar expansion 
for duu'{s,s^). The first two terms, which were explic- 
itly computed in lfT6ll for s, yield the "smooth part" 
{k — /^(5))/(27r), where ^ is the curvature (note that this 
can also be derived by applying the methods of [8] on 
the curved boundary corrections to Berry's random wave 
model, which are described in The third term is in- 

terpreted as a summation over possible paths from to s, 
where the nth summand includes integration over n inter- 
mediate bounce points. Approximating the integrals by the 
stationary phase method (closely following jlSn ). we get a 
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FIG. 1. The semi-classical and numerical length spectra (Fourier transform of d^{q)). (a): Absolute value, (b): Magnified view of 3 
prominent peaks. 



sum over classical paths (allowing only specular bounces). 
The oscillating part for duu' is given by 

where t enumerates classical orbits from to s, Lt is the 
length of the orbit, Pt is the Maslov index (number of con- 
jugate points plus twice the number of bounce points), ijj 
and are the angles between the orbit and the boundary 
at s and s' respectively, and p' = kco^{il)') is the clas- 
sical momentum at s' (hence, the oscillating part of duu' 
is O^s/k)). Taking derivatives of this expansion, we get 
an expression for duu' = dgds'duu' •> which has a similar 
form. Due to the rapid decay of the Fourier transformed 
convolution kernel p, the expansions (taken at 5 = s') 
of d^^^, and d^^^, converge, so the quotient required for 
substitution in © can be computed, and used to derive 
an expression for {b{s)). To compute df^{k), we multi- 
ply by d^{k) (for which we can use the Gutzwiller trace 
formula, again smoothed by p to ensure convergence), and 
integrate over s. The stationary phase condition added by 
this extra integration ensures that ij) = il^' , and the result 



includes summation over periodic classical orbits. Finally, 
we want to discard the spectral information and compute 
r}{n) rather than r]{k). Therefore, following the method 
described in f^], we substitute in the resulting expression 
for df^ (k), an expansion for /c(n), achieved by formally in- 
verting the Gutzwiller trace formula. This leads us to the 
trace formula for df^ (n), presented in ([T]). 

A stringent test of the theory above, is based on the fol- 
lowing argument. A "partial" trace formula which counts 
BI located on a prescribed part F C dVt of the boundary 
can be similarly derived by integrating the BI density over 
F alone: dj^r{k) = di){s; k)ds. Since the formula for 
di), much like (|7]), involves summation over orbits starting 
and ending at s, we conclude that only periodic orbits that 
have a bounce point in F will contribute to the sum in the 
resulting trace formula for d^r- By choosing a F which 
is bounded away from the bounce points of a specific orbit, 
we can effectively turn off the effect of that orbit. Similarly, 
we expect orbits that have some, but not all of their bounce 
points in the excluded regions dQ \ F, to have reduced am- 
plitude in the length spectrum. This result is demonstrated 
in FIG. El In FIG.Etb), F is plotted with a wide line, while 
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the excluded part dQ \ F is dotted. For demonstration pur- 
poses, the BI of ?/^i5o are shown on the boundary (total 
77150 = 24), and the points to be excluded from rjr are 
marked with empty circles (77r 150 = 16). Three orbits are 
shown, and the corresponding peaks in the length spectrum 
are also marked in FIG. Ha). Orbit 1 has both its bounce 
points in the excluded regions, so it completely disappears 
from the length spectrum corresponding to the partial count 
77r. Orbit 2, which has both of its bounce points in F is not 
effected by the exclusion, and orbit 3, which has only 1 out 
of 4 bounce points in F, is significantly inhibited, and drops 
below the noise level for the numerical case. This test and 
the general agreement between the semi-classical and the 
numerical length spectra give credence to the validity of 
the proposed trace formula. 
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